\input{../newcommands}

In \cite{Coleieee1997,Coleieee2002}, Cole introduced an implementation
of the source-free Maxwell's wave equations for narrow-band applications
based on non-standard finite-differences (NSFD). In \cite{Karkicap06},
Karkkainen \emph{et al.} adapted it for wideband applications. At
the Courant limit for the time step and for a given set of parameters,
the stencil proposed in \cite{Karkicap06} has no numerical dispersion
along the principal axes, provided that the cell size is the same
along each dimension (i.e. cubic cells in 3D). The ``Cole-Karkkainnen''
(or CK) solver uses the non-standard finite difference formulation
(based on extended stencils) of the Maxwell-Ampere equation and can be 
implemented as follows \cite{Vayjcp2011}:

\begin{subequations}
\begin{eqnarray}
D_{t}\mathbf{B} & = & -\nabla^{*}\times\mathbf{E}\label{Eq:Faraday}\\
D_{t}\mathbf{E} & = & \nabla\times\mathbf{B}-\mathbf{J}\label{Eq:Ampere}\\
\left[\nabla\cdot\mathbf{E}\right. & = & \left.\rho\right]\label{Eq:Gauss}\\
\left[\nabla^{*}\cdot\mathbf{B}\right. & = & \left.0\right]\label{Eq:divb}
\end{eqnarray}
\end{subequations}

Eq. \ref{Eq:Gauss} and \ref{Eq:divb} are not being solved explicitly
but verified via appropriate initial conditions and current deposition
procedure. The NSFD differential operators is given by $\nabla^{*}=D_{x}^{*}\mathbf{\hat{x}}+D_{y}^{*}\mathbf{\hat{y}}+D_{z}^{*}\mathbf{\hat{z}}$
where $D_{x}^{*}=\left(\alpha+\beta S_{x}^{1}+\xi S_{x}^{2}\right)D_{x}$
with $S_{x}^{1}G|_{i,j,k}^{n}=G|_{i,j+1,k}^{n}+G|_{i,j-1,k}^{n}+G|_{i,j,k+1}^{n}+G|_{i,j,k-1}^{n}$,
$S_{x}^{2}G|_{i,j,k}^{n}=G|_{i,j+1,k+1}^{n}+G|_{i,j-1,k+1}^{n}+G|_{i,j+1,k-1}^{n}+G|_{i,j-1,k-1}^{n}$.
$G$ is a sample vector component, while $\alpha$, $\beta$ and $\xi$
are constant scalars satisfying $\alpha+4\beta+4\xi=1$. As with
the FDTD algorithm, the quantities with half-integer are located between
the nodes (electric field components) or in the center of the cell
faces (magnetic field components). The operators along $y$ and $z$,
i.e. $D_{y}$, $D_{z}$, $D_{y}^{*}$, $D_{z}^{*}$, $S_{y}^{1}$,
$S_{z}^{1}$, $S_{y}^{2}$, and $S_{z}^{2}$, are obtained by circular
permutation of the indices.

Assuming cubic cells ($\Delta x=\Delta y=\Delta z$), the coefficients
given in \cite{Karkicap06} ($\alpha=7/12$, $\beta=1/12$ and $\xi=1/48$)
allow for the Courant condition to be at $\Delta t=\Delta x$, which
equates to having no numerical dispersion along the principal axes.
The algorithm reduces to the FDTD algorithm with $\alpha=1$ and $\beta=\xi=0$.
An extension to non-cubic cells is provided by Cowan, \emph{et al.}
in 3-D in \cite{CowanPRSTAB13} and was given by Pukhov in 2-D in
\cite{PukhovJPP99}. An alternative NSFDTD implementation that enables superluminous waves is also
given by Lehe {\it et al.} in \cite{LehePRSTAB13}. 

As mentioned above, a key feature of the algorithms based on NSFDTD
is that some implementations \cite{Karkicap06,CowanPRSTAB13} enable the time step $\Delta t=\Delta x$ along one or
more axes and no numerical dispersion along those axes. However, as
shown in \cite{Vayjcp2011}, an instability develops at the Nyquist
wavelength at (or very near) such a timestep. It is also shown in
the same paper that removing the Nyquist component in all the source
terms using a bilinear filter (see description of the filter below)
suppresses this instability.

